Glacier mass balance
We simplify the glacier mass balance to the following equation
\[ \Delta S = P-M = A_{imbal} \]
where the change of water storage (\(\Delta S\)) is equal to the precipitation (\(P\), assumed solid) minus the glacier melt (\(M\)). Typically melt excesses precipitation and we have negative imbalance ablation, indicating glacier storage loss. The glacier mass balance is calculated in annual time steps. We thereby refer to the hydrological year starting on October 1st of the previous year to take a full accumulation and ablation season into account.
Temperature index model
Hock (2003) describes several variants of the temperature index model for simulating glacier melt. The package implements the temperature index melt model described in (Hock 2003) in the function .
\[ M = \biggl\{ \begin{array}{l, l} 0, & T < T_{threshold} \\ f_{M} \cdot \left( T - T_{threshold} \right), & T >= T_{threshold} \end{array} \]
where \(M\) is the glacier melt in \(mm/d\), \(T\) is the daily average temperature in \(^{\circ} C\). The two parameters \(f_{M}\) and \(T_{\text{threshold}}\) refer to the melt factor and the threshold temperature above which glacier melt occurs and need to be calibrated. They have the units \(\frac{mm}{^{\circ} C \cdot d}\) and \(^{\circ} C\) respectively. Glacier melt is calculated in daily time steps.
Glacier volume development
As glaciers melt, their volume changes. This has to be taken into account for the long-term simulation of glacier discharge. To determine the initial glacier volume we use the area of the geometry of the RGI v6.0 data set and multiply it with the average thickness of the glacier. Glaciers larger than 1km2 we sub-divide into elevation bands of 100m altitude to account for elevation dependent temperature forcing. The large glaciers melt from the lowest elevation band to the highest elevation band whereby the glacier melt is subtracted from the glacier volume of lowest elevation band which is still glacierized. The small glaciers are not spatially discretized and thus they are melted homogeneously.
For annual time step \(t\), the evolution of the glacier volume is calculated as follows:
\[ A(t) = \text{glacierArea_RGIF}\bigl(V(t)\bigr) \]
\[ Q_{glacier}(t) = M(t) \cdot A(t) \]
\[ V(t+1) = V(t) + \Delta S = V(t) + \text{glacierImbalAbl}\bigl(M(t)\bigr) \]
Note that the imbalance ablation (or \(\Delta S\)) is negative for glacier mass loss and positive for glacier accumulation.
\(Q_{\text{glacier}}\) can be calibrated against glacier discharge derived from the Miles & Hugonnet data sets. The automated calibration is not yet included in the package.
The function is an empirical scaling function analogue the inverse of the scaling function derived by Erasov (1968) but based on the modern RGI v6.0 glacier geometries and the glacier thickness data set by Farinotti et al. (2019). The package implements volume-area and area-volume scaling functions based on both, Erashov and RGI-Farinotti data, allowing the estimation of glacier areas based on glacier volumes ( and ) and estimations of glacier volumes based on glacier areas ( and ).
The function is an empirical scaling function relating glacier imbalance ablation, i.e. the glacier storage loss glacier melt. It is derived from the glacier discharge data set by Miles et al. (2021) and the glacier thinning data set by Hugonnet et al. (2021).
Propagation of uncertainty
The initial glacier volume of each glacier is attributed an uncertainty of p/m 26%. This number is based on the average uncertainty of the glacier volume per RGI region reported in Farinotti et al. (2019) (Table 1). The error is the same for RGI regions 13 and 14. The uncertainty of the area of the RGI v6.0 glacier outlines is not known. We therefore assume that the error of the glacier volume stems to 50% from the estimation of the glacier thickness and to 50% from the glacier area. We further assume that the errors of the glacier area and glacier thickness are un-correlated and can thus estimate the uncertainties of the glacier area data and the glacier thickness data to be p/m 13% each.
For the non-linear relationships and , the standard deviation of the residuals of the fit was computed. The estimated error of the fit is assumed to be equal plus/minus twice the standard deviation of the residuals and yields 31% and 53% respectively. The residuals are not normally distributed and their actual distribution in unknown. Nevertheless, we apply the principles of linear error propagation for the sake of simplicity. We assume uncorrelated errors and therefore add the error of the function input to the error of the fit to obtain the error of the function output.
\[ \varepsilon_{V} = \varepsilon_{A} + \varepsilon_{\text{glacierVolume_RGIF}} = 0.26 + 0.31 = 0.57 \]
\[ \varepsilon_{A} = \varepsilon_{V} + \varepsilon_{\text{glacierArea_RGIF}} = 0.26 + 0.53 = 0.79 \]
We do not have an estimate for the error of the temperature index model. We assume a conservative relative error of 2, indicating that the estimated glacier melt is within a range of plus/minus 2 times it’s value.
The error of the imbalance ablation amounts to 73%. Assuming independent errors from the temperature index model for glacier melt and the scaling function between imbalance ablation and glacier melt, the error of the imbalance ablation amounts to approximately:
\[ \varepsilon_{Q_{\text{imb,melt}}} = \varepsilon_{Q_{\text{M}}} + \varepsilon_{\text{glacierImbalAbl}} = 2 + 0.73 = 2.73 \]
The errors stated above relate to initial estimates of the glacier areas, volumes, total ablation and imbalance ablation. The uncertainties are not propagated through each time step of the model as the errors depend on the errors of the previous time step, i.e. they are not uncorrelated over time. In any case, the uncertainty margins for any glacier melt modelling done with the presented method are large.
# Relative error estimates for the initial state in %
error_stats = tibble::tibble(
sV = 0.26, # Glacier volume
sA = 0.13, # Glacier area
sTh = 0.13, # Glacier thickness
sMelt = 2, # Glacier melt
sVolumeRGIF = 0.31, # Calculating glacier volume from glacier area
sAreaRGIF = 0.53, # Calculating glacier area from glacier volume
sImbal = 0.73, # Calculating glacier imbalance ablation from glacier melt
sQglgrowth = 1.6, # Calculating total ablation from glacier thinning
sQglmelt = 1.6) Demonstration
We demonstrate the above described method with the data from the Atbashy basin, a tributary of the Naryn river. As for the previous chapter you can download all required data from this link.
library(tmap)
library(timetk)
library(sf)
library(raster)
library(exactextractr)
library(tidyverse)
library(lubridate)
devtools::install_github("hydrosolutions/riversCentralAsia")
library(riversCentralAsia)
# Stochastic optimisation algorithm for model calibration
library(GA)
# Path to the data directory downloaded from the download link provided above.
# Here the data is extracted to a folder called atbashy_glacier_demo_data
data_path <- "../../atbashy_glacier_demo_data/"Spatial discretization
Glacier melt is calculated on elevation bands. The following code chunks show how to derive elevation bands on glaciers.
Glaciers smaller than 1 km2 are not sub-divided to limit the use of computational resources. The area threshold of 1 km2 is arbitrary but similar area thresholds have been used in the literature (e.g. Miles et al. (2021)).
library(sf)
library(raster)
library(exactextractr)
library(tidyverse)
devtools::install_github("hydrosolutions/riversCentralAsia")
library(riversCentralAsia)
data_path <- "../../atbashy_glacier_demo_data/"
# Define output file name
outfile <- paste0(data_path, "GIS/rgi_glaciers_atbaschy_el_bands.shp")
# Read in HRUs for the HBV models
hru <- st_read(paste0(data_path, "GIS/16076_HRU.shp"))
# Read in glacier geometries. Here we have already intersected the glacier
# outlines from RGI v6.0 with the outlines of the sub-basins of the Atbashy
# model in QGIS, so we know which glacier lies in which sub-basin. For the
# generation of the elevation bands, you can also use a simple RGI layer.
rgi <- st_read(paste0(data_path, "GIS/16076_Glaciers_per_subbasin.shp")) |>
dplyr::select(RGIId, Area, Zmed) # Area in km2, Zmed in masl (WGS84)
rgi_large_glaciers <- rgi |> filter(Area >= 1)
rgi_small_glaciers <- rgi |> filter(Area < 1)
dem <- raster(paste0(data_path, "GIS/16076_DEM.tif")) |>
mask(rgi_large_glaciers) |>
setMinMax()
rgi_large_glaciers <- rgi_large_glaciers |> st_transform(crs = 4326)
rgi_small_glaciers <- rgi_small_glaciers |> st_transform(crs = 4326)
## Elevation Band Zonation Parameters
band_interval <- 100 # in meters
holeSize_km2 <- 0.2 # cleaning holes smaller than that size
smoothFact <- .5 # level of band smoothing
demAggFact <- 1 # dem aggregation factor (carefully fine-tune this)
## Generate elevation bands for the large glaciers
hru_shp <- gen_basinElevationBands(
dem_PathN = dem,
demAggFact = demAggFact,
band_interval = band_interval,
holeSize_km2 = holeSize_km2,
smoothFact = smoothFact)
hru_shp_latlon <- hru_shp |>
st_transform(crs = 4326) |>
st_make_valid() # now latlon
## Intersect the elevation bands with the rgi layer
rgi_large_glaciers_el_bands <- st_intersection(rgi_large_glaciers, hru_shp_latlon)
## Concatenate with the small glaciers
rgi_elbands <- rgi_large_glaciers_el_bands |>
add_row(rgi_small_glaciers |>
mutate(layer = 1))
# Now we have a GIS layer with small glacier represented by 1 elevation band and
# large glaciers represented by several elevation bands
st_write(hru_out, outfile, append = FALSE)We can do some post-processing of the GIS layer to add an estimate of the glacier volume and the areas of the elevation bands.
# Read DEM
dem <- raster(paste0(data_path, "GIS/16076_DEM.tif"))
# Read in HRUs for the HBV models
hru <- st_read(paste0(data_path, "GIS/16076_HRU.shp"))## Reading layer `16076_HRU' from data source
## `C:\Users\marti\hydrosolutions Dropbox\Bea martibeatrice@gmail.com\marti\github\atbashy_glacier_demo_data\GIS\16076_HRU.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 40 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1105700 ymin: 4575416 xmax: 1197076 ymax: 4616616
## Projected CRS: WGS 84 / UTM zone 42N
# Read in glacier geometries. Here we have already intersected the glacier
# outlines from RGI v6.0 with the outlines of the sub-basins of the Atbashy
# model in QGIS, so we know which glacier lies in which sub-basin.
rgi <- st_read(paste0(data_path, "GIS/16076_Glaciers_per_subbasin.shp")) |>
# Area in km2, Zmed in masl (WGS84)
dplyr::select(RGIId, Area, Zmed, Zmax, Slope, Aspect) |>
distinct(RGIId, .keep_all = TRUE)## Reading layer `16076_Glaciers_per_subbasin' from data source
## `C:\Users\marti\hydrosolutions Dropbox\Bea martibeatrice@gmail.com\marti\github\atbashy_glacier_demo_data\GIS\16076_Glaciers_per_subbasin.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 135 features and 24 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 76.2421 ymin: 41.11981 xmax: 77.3168 ymax: 41.41837
## Geodetic CRS: WGS 84
# Read in the glacier elevation bands layer
rgi_elbands <- st_read(paste0(data_path,
"GIS/rgi_glaciers_atbaschy_el_bands.shp"),
quiet = TRUE) |>
st_transform(crs = crs(dem))
# Load a pre-processed raster file with glacier thickness in the Atbashy basin
# (see vignette [glaciers 01]{glaciers-01-intro.html} for details on the pre
# -processing)
glacier_thickness <- raster(
paste0(data_path,
"GLACIERS/Farinotti/pre-processed_glacier_thickness.tif"))
# Only keep the variables we need for this analysis
rgi_elbands <- rgi_elbands |>
dplyr::select(RGIId, Area, elvtn_b) |>
rename(Area_tot_glacier_km2 = Area) |>
mutate(ID = paste0(RGIId, "_", elvtn_b))
# Get mean elevation of each glacier/elevation band from DEM
rgi_elbands$z_masl <- exact_extract(dem, rgi_elbands, "mean", progress = FALSE)
# Update the glacier area within the basin boundaries
rgi_elbands$A_km2 <- as.numeric(st_area(rgi_elbands))*10^(-6)
glaciers <- unique(rgi_elbands$RGIId)
# Calculate the average glacier thickness for each elevation band.
rgi_elbands$thickness_m = exact_extract(glacier_thickness,
rgi_elbands, "mean", progress = FALSE)
rgi_elbands$thickness_m[is.na(rgi_elbands$thickness_m)] <- 0
tmap_mode("view")
tm_shape(dem, name = "DEM") +
tm_raster(n = 6,
palette = terrain.colors(6),
alpha = 0.8,
legend.show = TRUE,
title = "Elevation (masl)") +
tm_shape(rgi_elbands,
name = "RGI v6.0 in elevation bands") +
tm_polygons(col = "z_masl", lwd = 0.4,
title = "Mean glacier\nelevation [masl]",
palette = "Blues") +
tm_scale_bar(position = c("left", "bottom"))Observed glacier melt
Observed average glacier thinning rates by Hugonnet et al. (2021) are shown in the figure below. 30 out of 135 glaciers in the Atbashy basin show increase of glacier thickness between 2000 and 2020.
# Load the glacier thinning data set and filter it to the glaciers in the
# catchment of interest.
hugonnet <- read_csv(paste0(
data_path, "/GLACIERS/Hugonnet/dh_13_rgi60_pergla_rates.csv"),
show_col_types = FALSE)
hugonnet <- hugonnet |>
dplyr::filter(rgiid %in% rgi$RGIId) |>
tidyr::separate(period, c("start", "end"), sep = "_") |>
mutate(start = as_date(start, format = "%Y-%m-%d"),
end = as_date(end, format = "%Y-%m-%d"),
period = round(as.numeric(end - start, units = "days")/366))
glaciers_hugonnet <- rgi |>
left_join(hugonnet |> dplyr::select(rgiid, area, start, end, dhdt, err_dhdt,
dvoldt, err_dvoldt, dmdt, err_dmdt,
dmdtda, err_dmdtda, period),
by = c("RGIId" = "rgiid"))
# Add glacier discharge by Miles et al.
miles <- read_delim(paste0(data_path,
"GLACIERS/Miles/Miles2021_Glaciers_summarytable_20210721.csv"),
show_col_types = FALSE, delim = ",") |>
dplyr::filter(RGIID %in% rgi$RGIId & VALID == 1)
glaciers_hugonnet <- glaciers_hugonnet |>
left_join(miles |> dplyr::select(RGIID, totAbl, totAblsig, imbalAbl,
imbalAblsig),
by = c("RGIId" = "RGIID")) |>
mutate(Qgl_m3a = ifelse(is.na(totAbl), glacierDischarge_HM(dhdt), -totAbl))
observed_melt <- glaciers_hugonnet |>
st_drop_geometry() |>
dplyr::filter(period == 1) |>
dplyr::select(RGIId, start, dhdt, Qgl_m3a, Area) |>
mutate(Area_m2 = Area * 10^6) |>
dplyr::select(-Area) |>
distinct()
totQgl <- glaciers_hugonnet |>
st_drop_geometry() |>
dplyr::filter(period == 20) |>
mutate(Qgl_m3a = ifelse(Qgl_m3a < 0, 0, Qgl_m3a),
dVweqda_m3a = ifelse(dmdtda<0, -dmdtda, 0) * area) |>
summarise(Qgl_m3a = sum(Qgl_m3a),
dVweqda_m3a = sum(dVweqda_m3a)) |>
transmute(Qgl_mcma = Qgl_m3a * 10^(-6),
Qgl_m3s = Qgl_m3a / (60*60*24*365),
dVweqda_mcma = dVweqda_m3a * 10^(-6),
dVweqda_m3s = dVweqda_m3a / (60*60*24*365))
tmap_mode("view")
tm_shape(dem, name = "DEM") +
tm_raster(n = 6,
palette = terrain.colors(6),
alpha = 0.8,
legend.show = TRUE,
title = "Elevation (masl)") +
tm_shape(glaciers_hugonnet |> dplyr::filter(period == 20) |>
mutate(Qgl_mcma = Qgl_m3a * 10^(-6)),
name = "Glacier discharge\n[10^6 m3/a]") +
tm_polygons(col = "Qgl_mcma", lwd = 0.4,
title = "Qgl (10^6 m3/a)") +
tm_scale_bar(position = c("right", "bottom"))The average annual contribution of glacier melt to river runoff is 1.9198103 m3/s. The annual norm discharge at Atbashy gauging station is 16.6 m3/s which yields an average annual glacier contribution of about 11.5651221 % which is reasonable given the high median altitude of the basin of 3470 masl.
Forcing
There will be a separate Vignette to demonstrate how to prepare the forcing. For now, please refer to the methodology described in the CAHAM book. We provide you with the pre-processed forcing data and give you just a brief overview over the data source. We have sub-divided the glaciers into elevation bands to be able to reproduce the elevation-dependent temperature forcing on the glacier.
Historical temperatures
As meteorological data for high elevations in Central Asia is very scarce we use the CHELSA v2.1 data set (Karger et al. 2017). This is a global data set of forcing data for hydrolgical models based on ERA5 but corrected for precipitation biases in mountainous regions and especially suited for high elevations. The daily CHELSA forcing has been cut to the Central Region by the originator of the data set, D. Karger, WSL and extracted to the hydrological response units of the glaciers in the Atbashy basin by T. Siegfried.
hist_obs <- readRDS(file = paste0(data_path, "CLIMATE/hist_obs_glacier_tas.rds"))
# Plot the temperature time series for a given glacier/elevation band
hist_obs[, c(1, 2, 5, 10)] |>
pivot_longer(-date, names_to = "layer", values_to = "Temp_degC") |>
mutate(layer = as.numeric(gsub(pattern = "_", "",
str_extract(layer, "_\\d$")))) |>
plot_time_series(
.date_var = date,
.value = Temp_degC,
.color_var = layer,
.color_lab = "El. band",
.plotly_slider = TRUE,
.smooth = FALSE,
.line_size = 0.5
)Daily temperature time series of the lowest, a middle-elevation and the highest elevation band of glacier RGI60-13.08930 in the Atbashy basin. Data source: CHELSA v2.1.
Please note that, as expected, the daily temperatures are higher in lower elevation bands (lower elevation band numbers) than in higher elevation bands.
Future temperatures
Future temperature development per glacier or elevation band is extracted from the 3 CMIP6 GCM models with highest priorities for the region (Krasting et al. 2018; Yukimoto et al. 2019; Tang et al. 2019) downloaded from COPERNICUS. We take 4 socioeconomic scenarios (SSP) into account, covering 4 different emission scenarios (RCP): SSP 1 for RCP 2.6, short ssp126; SSP 2 for RCP 4.5, short ssp245; SSP 3 for RCP 7.0, short ssp370; SSP 5 for RCP 8.5, short ssp585. The temperatures of the climate models are bias corrected using the CHELSA data and a quantile mapping method (Gudmundsson et al. 2012; Gudmundsson 2016). The figure below shows the development of the temperature in the lowest elevation band of glacier RGI60-13.08930 as an example.
fut_sim <- readRDS(file = paste0(data_path, "CLIMATE/fut_sim_glacier_tas_qmapped.RDS"))
# Plot the temperature time series for a given glacier/elevation band
glacier <- "RGI60-13.08930_1"
# Extract the temperature for the selected glaciers for all GCMs and SSPs
scenarios <- names(fut_sim)
fut_temp <- NULL
for (scenario in scenarios) {
fut_temp <- rbind(fut_temp,
fut_sim[[scenario]] |>
dplyr::select(Date, all_of(glacier)) |>
mutate(Scenario = scenario))
}
fut_temp <- fut_temp |>
mutate(Hyear = hyear(Date)) |>
group_by(Hyear, Scenario) |>
summarise(Date = first(Date),
Temp = mean(get(glacier))) |>
separate(Scenario, into = c("GCM", "SSP"), sep = "_") |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear),
GCM != "IPSL-CM6A-LR")
# Plot annual data for better readability
ggplot() +
geom_line(data = hist_obs |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear) |>
summarise(date = first(date),
Temp = mean(get(glacier))) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)),
aes(date, Temp)) +
geom_line(data = fut_temp, aes(Date, Temp, colour = GCM,
linetype = SSP)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "T [deg C]", title = glacier) +
theme_bw()
Annual historical (CHELSA) and future (CMIP6, Copernicus) temperature for the lowest elevation band of glacier RGI60-13.08930.
Calculate daily glacier melt
The code snippet below illustrates how to calculate the daily glacier melt in a river catchment with default parameters for the temperature index model. The figure shows daily melt rates for the different elevation bands of one of the larger glaciers (RGI60-13.08930). The highest melt occurs at low elevation bands. Glacier melt only occurs during the summer months when the average daily temperature is above the threshold temperature for glacier melt.
# Default parameters of the temperature index model for glacier melt.
MF_small = 2 # Melt coefficient for small glaciers (smaller than 1km2)
MF_large = 0.1 # Melt coefficient for large glaciers (larger or equal than 1km2)
threshold_temperature = 0 # Default temperature threshold for glacier melt
Area <- rgi_elbands |>
st_drop_geometry() |>
dplyr::select(ID, A_km2) |>
pivot_wider(names_from = ID, values_from = A_km2)
# Assign different melt factors to large and small glaciers.
MF <- Area |>
mutate(across(everything(), ~MF_large),
across(ends_with("_1"), ~MF_small))
melt <- glacierMelt_TI(temperature = hist_obs |> dplyr::select(-date),
MF = MF,
threshold_temperature = threshold_temperature)
melt <- as_tibble(melt) |>
mutate(date = hist_obs$date) |>
relocate(date, .before = where(is.numeric))
ggplot(melt[, 1:9] |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
left_join(rgi_elbands |>
st_drop_geometry() |>
dplyr::select(ID, z_masl), by = "ID") |>
mutate("Mean\nelevation\n[masl]" = as.factor(round(z_masl)))) +
geom_line(aes(date, Melt, colour = `Mean\nelevation\n[masl]`)) +
scale_colour_viridis_d() +
labs(x = "Date", y = "Melt [mm/d]") +
theme_bw()
Daily glacier melt per elevation band.
Compare to measured glacier melt
Aggregate the daily glacier melt to annual data and compare it to the observed glacier melt.
melt_a_eb <- melt |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear, ID) |>
summarise(date = first(date),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt))|>
ungroup()
melt_a <- melt_a_eb |>
separate(ID, into = c("RGIId", "elB"), sep = "_") |>
group_by(Hyear, RGIId) |>
summarise(date = as_date(first(date)),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt)) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)) |>
left_join(rgi |> dplyr::select(RGIId, Area) |> st_drop_geometry(),
by = "RGIId") |>
# Assume constant glacier area over the comparison period.
mutate(Qgl_m3a = Melt*10^(-3) * Area*10^6,
.lb_Qgl_m3a = Qgl_m3a * (1-sqrt(error_stats$sMelt^2+error_stats$sA^2)),
.lb_Qgl_m3a = ifelse(.lb_Qgl_m3a < 0, 0, .lb_Qgl_m3a),
.ub_Qgl_m3a = Qgl_m3a * (1+sqrt(error_stats$sMelt^2+error_stats$sA^2)))
melt_obs_a_all <- observed_melt |>
mutate(.lb_Qgl_m3a = ifelse(
dhdt > 0,
ifelse(Qgl_m3a*(1-error_stats$sQglgrowth)<0, 0,
Qgl_m3a*(1-error_stats$sQglgrowth)),
ifelse(Qgl_m3a*(1-error_stats$sQglmelt)<0, 0,
Qgl_m3a*(1-error_stats$sQglmelt))),
.ub_Qgl_m3a = ifelse(dhdt > 0,
Qgl_m3a*(1+error_stats$sQglgrowth),
Qgl_m3a*(1+error_stats$sQglmelt)))
melt_obs_a <- melt_obs_a_all |>
dplyr::filter(RGIId %in% glaciers[c(6, 7)])
ggplot() +
geom_ribbon(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[c(6,7)]),
aes(date, Qgl_m3a, ymin = .lb_Qgl_m3a, ymax = .ub_Qgl_m3a,
fill = RGIId), colour = NA, alpha = 0.1) +
geom_ribbon(data = melt_obs_a,
aes(start, Qgl_m3a,
ymin = .lb_Qgl_m3a,
ymax = .ub_Qgl_m3a,
colour = RGIId, linetype = "obs", fill = RGIId),
size = 0.2, alpha = 0.1) +
geom_line(data = melt_a |>
dplyr::filter(RGIId %in% glaciers[c(6,7)]),
aes(date, Qgl_m3a, colour = RGIId, linetype = "sim")) +
geom_line(data = melt_obs_a, aes(start, Qgl_m3a, colour = RGIId,
linetype = "obs")) +
labs(x = "Date", y = "Melt [m weq/a]") +
scale_linetype_manual(name = "Source",
values = c("sim" = 1, "obs" = 2)) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ggtitle(paste0("MF: ", MF, "Tth: ", threshold_temperature)) +
theme_bw()
Daily glacier melt per elevation band.
The temperature index model can be calibrated with the specific glacier volume change provided by Hugonnet et al. (2021) (see \(\ref{fig:glaciermeltcomparison}\)).
Calibration of the temperature index model
The 2 parameters of the temperature index model can be calibrated using observed glacier melt rates. As mentioned above, the observed glacier melt rates are simulated total ablation by Miles et al. (2021) where available or estimated total ablation using an empirical relationship between total glacier ablation by Miles et al. (2021) and glacier thinning rates by Hugonnet et al. (2021).
The calibration is done using an R implementation of a genetic
algorithm Scrucca (2017). The optimization
is slow but it is necessary to use a global optimizer
(nlminb of the stats package for example will not converge
for all glaciers).
# Calibrating many glaciers takes a lot of time.
# Reduce the amount of time required and only calibrate the glaciers that are
# larger or equal than 0.2 km2. That represents 83% of the glacier volume in the
# Atbashy basin but reduces the number of glaciers to be calibrated by 65%.
# Or, in the cas of an Area threshold of 0.5 km2, the glacier volume covered is
# 59% and the glacier number reduction is 92%
rgi <- rgi |> mutate(Vol_km3 = glacierVolume_RGIF(Area))
rgi_cal <- rgi |>
dplyr::filter(Area >= 0.1)
temperature_cal <- hist_obs |>
mutate(date = as_date(date)) |>
dplyr::select(c(date, starts_with(rgi_cal$RGIId)))
observed_cal <- observed_melt |>
rename(RGIID = RGIId,
totAbl = Qgl_m3a,
year = start) |>
mutate(year = year(year)) |>
dplyr::filter(RGIID %in% rgi_cal$RGIId)
glaciers_cal <- rgi_cal$RGIId
DDFfit <- NULL
for (id in c(1:length(glaciers_cal))) {
cat("Fitting temperature index model for glacier", glaciers_cal[id], "id =",
id, "\n")
tempGA <- ga(type = "real-valued",
fitness = glacierDischargeRMSE,
temperature = temperature_cal,
observed = observed_cal,
index = id,
lower = c(0.1, 0), upper = c(50, 0),
popSize = 500, maxiter = 100,
parallel = (parallel::detectCores()-6))
cat(" Result: RMSE = ", -tempGA@fitnessValue, "MF =", tempGA@solution[1],
"Tth =", tempGA@solution[2], "\n")
DDFfit <- rbind(DDFfit, tibble(MF = tempGA@solution[1],
Tth = tempGA@solution[2],
rsme = tempGA@fitnessValue) |>
mutate(ID = glaciers_cal[id]))
save(DDFfit,
file = paste0(data_path, "GLACIERS/TIparameters.RData")
}How to determine the MF parameter for the non-calibrated glaciers? Using the average MF from the calibrated glaciers is one option. The best results (for the present case study) are obtained by modulating the MF according to glacier slope.
Now we re-calculate the glacier melt using the calibrated data set.
# Option to calculate only for a selection of glaciers (modify indices)
indices <- c(1:length(glaciers))
DDFfit |> dplyr::filter(RGIId %in% glaciers[indices])## # A tibble: 179 x 9
## MF Tth rsme RGIId ID Slope .fitted .resid layer
## <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_1 23.7 4.84 -3.93 1
## 2 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_2 23.7 4.84 -3.93 2
## 3 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_3 23.7 4.84 -3.93 3
## 4 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_4 23.7 4.84 -3.93 4
## 5 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_5 23.7 4.84 -3.93 5
## 6 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_6 23.7 4.84 -3.93 6
## 7 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_7 23.7 4.84 -3.93 7
## 8 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_8 23.7 4.84 -3.93 8
## 9 0.916 0 -115. RGI60-13.06145 RGI60-13.06145_9 23.7 4.84 -3.93 9
## 10 0.547 0 -72.6 RGI60-13.06157 RGI60-13.06157_1 22.8 4.77 -4.23 1
## # ... with 169 more rows
melt_cal <- glacierMelt_TI(temperature = hist_obs |>
dplyr::select(-date),
MF = MFvec,
threshold_temperature = Tthvec) |>
mutate(date = hist_obs$date) |>
relocate(date, .before = where(is.numeric))
melt_a_eb <- melt_cal |>
pivot_longer(-date, names_to = "ID", values_to = "Melt") |>
mutate(Hyear = hyear(date)) |>
group_by(Hyear, ID) |>
summarise(date = first(date),
Melt = sum(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt))|>
ungroup()
melt_a <- melt_a_eb |>
separate(ID, into = c("RGIId", "elB"), sep = "_") |>
group_by(Hyear, RGIId) |>
summarise(date = as_date(first(date)),
Melt = mean(Melt),
.lb_Melt = ifelse(Melt*(1-error_stats$sMelt)<0, 0,
Melt*(1-error_stats$sMelt)),
.ub_Melt = Melt*(1+error_stats$sMelt)) |>
ungroup() |>
dplyr::filter(Hyear > min(Hyear) & Hyear < max(Hyear)) |>
left_join(rgi |> dplyr::select(RGIId, Area) |> st_drop_geometry(),
by = "RGIId") |>
# Assume constant glacier area over the comparison period.
mutate(Qgl_m3a = Melt*10^(-3) * Area*10^6,
.lb_Qgl_m3a = Qgl_m3a * (1-sqrt(error_stats$sMelt^2+error_stats$sA^2)),
.lb_Qgl_m3a = ifelse(.lb_Qgl_m3a < 0, 0, .lb_Qgl_m3a),
.ub_Qgl_m3a = Qgl_m3a * (1+sqrt(error_stats$sMelt^2+error_stats$sA^2)))
melt_a_basin <- melt_a |>
group_by(Hyear) |>
summarise(Qgl_m3a = sum(Qgl_m3a, na.rm = TRUE),
.lb_Qgl_m3a = sum(.lb_Qgl_m3a, na.rm = TRUE),
.ub_Qgl_m3a = sum(.ub_Qgl_m3a, na.rm = TRUE),
date = first(date),
Area = sum(Area, na.rm = TRUE),
Melt = mean(Melt, na.rm = TRUE)) |>
ungroup()
melt_obs_a_basin <- melt_obs_a_all |>
mutate(Hyear = hyear(start)) |>
group_by(Hyear) |>
summarise(Qgl_m3a = sum(Qgl_m3a, na.rm = TRUE),
.lb_Qgl_m3a = sum(.lb_Qgl_m3a, na.rm = TRUE),
.ub_Qgl_m3a = sum(.ub_Qgl_m3a, na.rm = TRUE),
date = first(start),
Area = sum(Area_m2, na.rm = TRUE)) |>
ungroup()
ggplot() +
geom_ribbon(data = melt_a_basin,
aes(date, Qgl_m3a, ymin = .lb_Qgl_m3a, ymax = .ub_Qgl_m3a),
colour = NA, alpha = 0.1) +
geom_ribbon(data = melt_obs_a_basin,
aes(date, Qgl_m3a,
ymin = .lb_Qgl_m3a,
ymax = .ub_Qgl_m3a, linetype = "obs"),
size = 0.2, alpha = 0.1) +
geom_line(data = melt_a_basin,
aes(date, Qgl_m3a, linetype = "sim")) +
geom_line(data = melt_obs_a_basin,
aes(date, Qgl_m3a, linetype = "obs")) +
labs(x = "Date", y = "Melt [m weq/a]") +
scale_linetype_manual(name = "Source",
values = c("sim" = 1, "obs" = 2)) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
theme_bw()
Glacier discharge in the Atbashy basin.
Aggregated to the basin level, the simulated glacier discharge fits the observed glacier discharge quite nicely.
comparison <- melt_obs_a_all |>
dplyr::filter(RGIId %in% glaciers[indices]) |>
mutate(Hyear = hyear(start)) |>
rename(Obs = Qgl_m3a,
.lobs = .lb_Qgl_m3a,
.uobs = .ub_Qgl_m3a) |>
left_join(melt_a |>
dplyr::filter(RGIId %in% glaciers[indices]) |>
rename(Sim = Qgl_m3a,
.lsim = .lb_Qgl_m3a,
.usim = .ub_Qgl_m3a),
by = c("Hyear", "RGIId"))
axmin <- min(cbind(comparison$Obs, comparison$Sim), na.rm = TRUE)
axmax <- max(cbind(comparison$Obs, comparison$Sim), na.rm = TRUE)
ggplot(comparison |> drop_na()) +
geom_abline(slope = 1, intercept = 0) +
geom_point(aes(Obs, Sim, colour = RGIId), size = 0.4, show.legend = FALSE) +
scale_colour_viridis_d() +
xlim(axmin, axmax) +
ylim(axmin, axmax) +
coord_fixed() +
theme_bw()
Observed versus simulated glacier discharge. 11 years of annual data for 131 glaciers in the basin.
However, when looking at calibration outcomes of individual glaciers, the picture does not look quite as well. For some glaciers, the optimization does not converge to produce good estimates of glacier discharge.
Glacier mass balance
The following figure shows the components of the glacier mass balance for a few glaciers in the Atbashy basin.
glacier_balance <- glacierBalance(melt_a_eb = melt_a_eb,
rgi_elbands = rgi_elbands)
glacier_balance <- glacier_balance |>
mutate(.lb = ifelse(Variable == "A_km2",
ifelse(Value*(1-error_stats$sA)>0,
Value*(1-error_stats$sA), 0),
ifelse(Variable == "V_km3",
ifelse(Value*(1-error_stats$sV)>0,
Value*(1-error_stats$sV), 0),
ifelse(Variable == "Q_m3a",
ifelse(Value*(1-error_stats$sQglmelt)>0,
Value*(1-error_stats$sQglmelt), 0),
ifelse(Variable == "Qimb_m3a",
Value*(1-error_stats$sImbal), NA)))),
.ub = ifelse(Variable == "A_km2",
Value*(1+error_stats$sA),
ifelse(Variable == "V_km3",
Value*(1+error_stats$sV),
ifelse(Variable == "Q_m3a",
Value*(1+error_stats$sQglmelt),
ifelse(Variable == "Qimb_m3a",
Value*(1+error_stats$sImbal), NA)))))
ggplot(glacier_balance |>
dplyr::filter(RGIId %in% glaciers[6:7],
Hyear > min(Hyear) & Hyear < max(Hyear),
Variable %in% c("A_km2", "V_km3", "Q_m3a", "Qimb_m3a"))) +
geom_ribbon(aes(Hyear, ymin = .lb, ymax = .ub, fill = RGIId),
alpha = 0.2, colour = NA) +
geom_line(aes(Hyear, Value, colour = RGIId)) +
facet_wrap("Variable", scales = "free_y") +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
theme_bw()
Glacier development and discharge.
We now have glacier discharge (Q_m3s) and the unsustainable contribution to glacier discharge, the imbalance ablation (Qimb_m3a) which is negative for glacier loss and positive for growing glaciers.
Future glacier melt
scenarios <- names(fut_sim)
fut_melt_d <- NULL
fut_melt_a_eb <- NULL
fut_glacier_balance <- NULL
for (scenario in scenarios) {
temp_melt_d <- glacierMelt_TI(temperature = fut_sim[[scenario]] |>
dplyr::select(-Date),
MF = MFvec,
threshold_temperature = Tthvec) |>
mutate(Scenario = scenario,
date = fut_sim[[scenario]]$Date) |>
relocate(date, .before = where(is.numeric)) |>
pivot_longer(-c(date, Scenario), names_to = "ID",
values_to = "Melt")
temp_melt <- temp_melt_d |>
mutate(Hyear = hyear(date)) |>
# Sum melt per year
group_by(Hyear, Scenario, ID) |>
summarise(date = first(date),
Melt = sum(Melt)) |>
ungroup()
fut_melt_d <- rbind(fut_melt_d, temp_melt_d)
fut_melt_a_eb <- rbind(fut_melt_a_eb, temp_melt)
fut_glacier_balance <- rbind(
fut_glacier_balance,
glacierBalance(melt_a_eb = temp_melt,
rgi_elbands = rgi_elbands) |>
mutate(Scenario = scenario))
}
fut_melt_amean <- fut_melt_a_eb |>
# Mean melt rate per glacier
# mutate(RGIId = gsub("\\_\\{d}+$", "", ID)) |>
# group_by(Hyear, Scenario, RGIId) |>
# summarise(Melt = mean(Melt)) |>
# Mean melt rate in the catchment
group_by(Hyear, Scenario) |>
summarise(date = as_date(first(date)),
Melt_mma = mean(Melt)) |>
ungroup() |>
separate(Scenario, into = c("GCM", "SSP"), sep = "_") |>
mutate(SSP = factor(SSP, levels = c("ssp585", "ssp370", "ssp245", "ssp126"))) |>
# dplyr::filter(!str_detect(GCM, "UK")) |>
dplyr::filter(!str_detect(GCM, "IPSL"),
Hyear < max(Hyear)) |>
group_by(Hyear, SSP) |>
mutate(lb = min(Melt_mma),
ub = max(Melt_mma),
mean = mean(Melt_mma)) |>
ungroup()
melt_amean <- melt_a |>
dplyr::select(Hyear, date, Melt) |>
group_by(Hyear) |>
dplyr::summarise(date = first(date),
Melt_mma = mean(Melt)) |>
ungroup()
ggplot(fut_melt_amean) +
geom_ribbon(aes(date, ymin = lb, ymax = ub, fill = SSP), alpha = 0.2) +
geom_line(aes(date, mean, colour = SSP)) +
geom_line(data = melt_amean, aes(date, Melt_mma)) +
# geom_smooth(aes(date, mean, colour = SSP), method = "loess", se = FALSE) +
# geom_smooth(data = melt_amean, aes(date, Melt_mma), method = "loess",
# se = FALSE) +
scale_colour_viridis_d() +
scale_fill_viridis_d() +
ylab("Melt [mm/a]") +
xlab("Year") +
theme_bw()
Future specific glacier melt rates expected in the Atbashy basin.
We see an increase of the expected glacier melt rates. In the two lesser emission scenarios ssp126 and ssp245, glacier melt rates stabilise around mit century and in the last quarter of the century repsectively on a higher level than in recent years while under the higher emission scenarios ssp370 and ssp585 glacier melt continues to rise.
Uncertainties
The shaded area in the above figure indicates the spread of the GCM ensemble used. The width of the uncertainty band of the model method is at least twice the width of the above uncertainty bands, such that it is near impossible to differentiate between the impact of the different socio-economic scenarios. Dramatic glacier volume reduction cannot be avoided even with the least emissions scenario ssp126.
Next steps
Now we have glacier discharge and glacier imbalance ablation. How do we include these in our modelling workflow? Glacier vignette 4 explains how glaciers are treated in the hydraulic/hydrologic model RS Minerve and details a workflow to implement glacier discharge in a model in RS Minerve.
Open Global Glacier Model (OGGM)
The Open Global Glacier Model is a community project which provides open-source glacier modelling tools in python and educational resources. The modeling framework includes comprehensive tutorials and data-pre and post-processing tools. It can be used as an alternative to the glacier modelling work flow presented here.